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This  study  focused  on  the  target  state  estimator, 
implemented  as  an  extended  Kalman  filter,  on  the  AFTI/F-16. 
Several  possible  reasons  for  the  poor  performance  it  exhibits 
are  investigated;  the  major  reason  for  this  poor  performance 
is  shown  to  lie  chiefly  with  the  conventional  covariance 
update  techniques  it  uses.  Others  have  shown  that  the 
recursion  equations  which  use  conventional  update  techniques 
are  numerically  unstable;  because  alternatives  exist,  the 
conventional  updates  should  never  be  used. 
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Abstract 

The  purpose  o£  this  study  was  to  investigate  target 
state  estimation  techniques  for  the  air-to-air  mode  of  the 
AFTl/F-16  automated  maneuvering  and  attack  system.  The 
target  state  estimator  (TSE)  previously  developed  would  not 
perform  to  specif ications;  possible  reasons  for  this  poor 
performance  are  presented  as  well  as  suggestions  to  upgrade 
the  performance. 

The  TSE  exists  as  extended  Kalman  filter  equations  in  a 
digital  computer.  The  previously  developed  Kalman  filter 
equations  used  conventional  covariance  update  techniques  and 
a  Gauss-Markov  system  dynamics  model  which  expressed  the 
states  in  an  inertial  reference  frame.  Measurements  were 
performed  in  the  line-o£-sight  (LOS)  frame,  but  the 
covariance  matrix  was  not  rotated  into  the  LOS  frame  during 
update.  This  study  focused  on  three  areas:  (1)  Determine  if 
the  Gauss-Markov  dynamics  model  was  adequate  for  the  tracking 
accuracies  specified.  (2)  Determine  if  a  rotation  had  to  be 
performed  to  account  for  the  states  being  expressed  in  one 
frame  while  the  measurements  were  physically  made  in  another. 
(3)  Determine  what  effect  the  conventional  covariance 
updates,  coupled  with  the  short  (16-bit)  wordlength  of  the 
TSE  computers,  has  on  the  stability  of  the  Kalman  filter. 

Two  filter  dynamics  models  were  designed,  tested,  and 
compared.  The  first  model  used  complex  equations  and  closely 
modeled  an  air-to-air  engagement.  Most  of  the  complexity  of 


the  model  was  maintained  in  its  implementation,  and  it  was 
used  as  a  baseline  model.  The  second  filter  used  a  Gauss- 
Markov  dynamics  model  and  made  several  assumptions  to 
simplify  computations. 

Analysis  of  filter  performance  revealed  that  the  Gauss- 
Markov  filter  dynamics  model  was,  indeed,  an  adequate  model. 
Also,  the  covariance  matrix  does  not  have  to  be  rotated  into 
the  LOS  frame  if  the  measurements  are  redefined.  The  second 
filter  was  then  implemented  using  U-D  covariance 
factorization  algorithms,  but  the  time  propagation  routines 
used  were  apparently  flawed.  However,  the  poor  performance 
of  the  TSE  is  no  doubt  caused  by  the  conventional  Kalman 
filter  recursions,  as  they  are  inherently  unstable. 
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IMPLEMENTATION  OF  A  TARGET  STATE  ESTIMATOR  FOR 
THE  AIR-TO-AIR  ATTACK  MODE  OF  THE  AFTI/F-16 


Introduction 


Background 

The  Flight  Dynamics  Laboratory  ( AFWAL/FI GX )  is 
evaluating  software  for  the  Automated  Maneuvering  and  Attack 
System  (AMAS)  for  the  Advanced  Fighter  Technology  Integration 
AFTI/F-16  aircraft.  The  AMAS  uses  integrated  fire/flight 
control  ( IFFC )  techniques  to  track  a  designated  target, 
calculate  a  flight  condition  from  which  a  released  weapon  has 
a  high  kill  probability,  maneuver  the  aircraft  to  that  flight 
condition,  and  release  the  weapon  [21. 

The  AFTI/F-16  is  capable  of  attacking  ground  as  well  as 
airborne  targets;  the  AMAS  should  automate  an  attack  against 
either.  Because  target  maneuverability  has  a  large  affect  on 
the  tracking  system's  capabilities,  the  AMAS  is  subdivided 
into  three  modes:  a  stationary  ground  target  mode,  for  use 
against  structures  or  other  non-mobile  targets;  a  mobile 
ground  target  mode,  for  use  against  moving  vehicles;  and  an 
airborne  target  mode,  capable  for  use  in  an  air-to-air 
engagement.  The  ground  modes  of  the  AMAS  are  currently 
implemented  in  the  AFTI/F-16.  However,  the  airborne  mode  is 
still  under  evaluation,  partly  because  the  AMAS  tracking 
function  does  not  always  perform  with  the  precision  required 
to  obtain  a  gun  firing  solution  [2,10]. 


The  AMAS  tracking  function  consists  of  sensors,  which 
detect  the  relative  position  and  velocity  of  the  target,  and 
a  target  state  estimator  (TSE),  which  estimates  current 
target  states.  The  sensors  are  a  forward-looking  infrared 
( FLIR )  imaging  system  and  the  inertial  navigation  unit  (INU) 
for  angle  measurements,  a  laser  range  finder,  and  the  APG-66 
radar  for  range  rate  (Doppler)  information.  The  TSE  consists 
of  an  extended  Kalman  filter  ( EKF )  algorithm  implemented  in  a 
digital  computer.  The  EKF  calculates  optimal  (as  described 
below)  estimates  of  current  target  states  [101. 

The  Fire  Control  Computer  (FCC)  performs  the  AMAS 
function  of  calculating  the  desired  flight  profile.  A 
discussion  of  the  algorithms  used  by  the  FCC  to  compute  the 
desired  flight  profile  is  beyond  the  scope  of  this  thesis. 
However,  the  knowledge  that  the  FCC  must  be  able  to  predict 
the  target  position  and  velocity  at  a  specific  time  in  the 
future  implies  that  the  target  states  must  Include  position, 
velocity,  and  acceleration.  To  ensure  that  body  rates  do  not 
influence  the  target  states,  the  target  states  are  expressed 
in  an  inertial  Cartesian  coordinate  system.  Therefore,  the 
states  estimated  by  the  TSE  are  relative  (target  with  respect 
to  attacker)  position,  relative  velocity,  and  total  target 
acceleration  expressed  in  an  inertial  Cartesian  coordinate 
system  (101. 

The  TSE  uses  an  EKF  to  combine  the  measurements 
available  from  the  sensors  with  a  target  dynamics  model  to 


estimates  are  the  best  estimates  possible  if  the  model  used 
in  the  filter  adequately  represents  the  real  world  system  and 
the  variances  used  in  the  filter  are  truly  representative  of 
the  variances  in  the  real  world  system  [4:41. 

Theoretically,  linear  Kalman  filters  provide  the  best 
state  estimates  possible  if  the  conditions  described  above 
are  met;  however,  the  algorithms  for  their  implementation  may 
not  be  numerically  stable.  This  numerical  instability  is 
also  apparent  for  some  implementations  of  EKFs .  In  fact,  the 
finite  wordlength  of  digital  computers  may  make  the 
recursion  equations  of  the  Kalman  filter  diverge,  as  they  do 
for  certain  conditions  of  the  currently  implemented  AFTI/F-16 
TSE  [2],  The  divergent  effects  of  finite  wordlength  may  be 
reduced  by  increasing  the  wordlength  of  the  computer  or  by 
modifying  the  recursion  equations  which  comprise  the  Kalman 
filter  [4:368],  Computer  wordlengths  are  difficult  to 
change,  whereas  the  recursion  equations  in  software  are 
readily  modified.  Possible  modifications  to  the  recursion 
algorithms  are  discussed  further  in  Chapter  II. 


1.2.  Er.oJblem  Statement 

The  EKF  currently  implemented  in  the  air-to-air  mode  of 
the  AFTI/F-16  TSE  diverges  for  certain  attack  geometries, 
causing  the  AMAS  to  break  off  the  attack  [21.  The  specific 
geometries  where  the  divergence  occurs  were  not  available  to 
the  researcher,  but  two  possible  causes  for  this  divergence 


are  addressed  by  this  research.  First,  the  current  EKF  uses 


standard  covariance  matrix  update  equations,  which  have  been 
shown  to  be  numerically  unstable  [1:338].  A  modification  to 
the  EKF  recursion  algorithms  is  proposed  by  this  thesis  to 
eliminate,  or  at  least  significantly  reduce,  this  problem. 
Second,  the  measurement  covariances  are  known  only  in  the 
line-of-sight  (LOS)  frame  while  the  target  states  are  known 
in  the  Inertial  frame.  Consequently,  a  change  of  reference 
frames  must  be  accomplished  before  the  states  and  covariances 
are  updated,  but  this  is  not  done  in  the  current  TSE  [10]. 

The  rotation  matrix  to  accomplish  this  change  of  reference 
frame  and  its  use  are  discussed  further  in  Chapter  II. 

1 .3.  Scope  and  Limitations 

This  thesis  develops  an  EKF  based  upon  the  model  used  in 
the  current  AFTI/F-16  TSE  and  compares  it  to  an  EKF  based 
upon  a  more  complex  model.  The  filters  developed  herein  use 
numerically  stable  update  algorithms,  reducing  or  eliminating 
divergence  due  to  the  digital  computer's  finite  wordlength. 
Several  simplifying  assumptions  are  made  which  significantly 
reduce  the  complexity  of  the  development,  as  discussed  below. 

The  first  assumption  is  that  the  measurement  coordinate 
frame  (the  LOS  frame)  is  assumed  to  be  inertially  space- 
stabilized  during  the  measurement.  While  this  is  not 
generally  true  of  AFTI/F-16  hardware  [10]  (the  sensor  head 
rotates  with  the  relative  target  position,  as  explained  in 
Chapter  II),  any  rotation  is  assumed  to  have  second  or  higher 
order  effects  on  TSE  performance.  This  assumption  eliminates 
the  requirement  to  calculate  Coriolis  acceleration  effects  on 


the  measurement;  instead,  errors  introduced  by  Coriolis 
accelerations  are  assumed  to  be  negligible  compared  to  the 
dynamics  noise  models  of  the  EKFs  . 

The  next  assumption  is  that  the  attacker's  position, 
velocity,  acceleration,  and  attitude  are  available  from  an 
inertial  navigation  unit  (INU)  without  error.  This 
assumption  is  made  because  current  INUs  have  errors  much 
smaller  than  the  errors  expected  from  the  TSE  [8:4].  Also, 
because  any  INU  errors  are  also  present  in  fire  control 
computer  calculations,  the  overall  effect  of  INU  errors  on 
the  AMAS  solution  is  reduced  [101. 

Another  assumption  is  that  the  parallax  errors  which 
occur  because  of  the  physical  separation  of  the  STS  and  the 
INU  are  negligible  compared  to  other  system  errors.  This  is 
justifiable  because  the  separations  are  small  (less  than  10 
feet)  : 10 ] , 

A  further  assumption  is  that  the  earth-fixed  (inertial) 
reference  frame  has  its  positive  axes  oriented  in  the  local 
north,  east,  and  down  (towards  the  center  of  the  earth) 
directions.  The  origin  of  this  coordinate  frame  is  fixed  at 
the  center  of  gravity  (eg)  of  the  attacker  at  the  instant 
target  designation  takes  place  (when  the  TSE  is  at  time 
zero )  . 

The  final  major  assumption  is  that  the  attacker  is 
benign;  l.e.,  the  attacking  aircraft  does  not  maneuver.  Thl 
assumption  is  made  to  simplify  the  implementation  of  testing 
routines.  Formulating  a  complete  model  of  the  AFTI/F-16 


control  system  and  aerodynamics  requires  more  time  than  is 
available  to  complete  this  thesis,  and  the  TSE  relies  on 
other  aircraft  systems  only  for  input  of  data.  Further,  the 
TSE  does  not  directly  control  any  aircraft  function. 

LJL.  Equipment  Used 

Each  of  the  Kalman  filter  implementations  investigated 
in  this  thesis  is  evaluated  by  means  of  software  written  for 
this  purpose,  the  Simulation  for  Optimal  Filter  Evaluation 
( SOFE )  and  associated  plotting  post-processor  (SOFEPL) 
routines  [6,7].  Each  of  these  tools  are  available  on  the 
AFIT/ENG  ELXSI  computer.  A  brief  discussion  of  the 
evaluation  process  follows;  a  more  detailed  description  of 
SOFE  and  SOFEPL  is  presented  in  Chapters  1 1  to  v. 

SOFE  is  an  EKF  evaluation  tool  which  was  written  to 
test  the  effects  of  varying  the  parameters  or  equations  which 
comprise  a  given  Kalman  filter.  Written  in  FORTRAN,  SOFE 
consists  of  a  set  of  fixed  subroutines  that  exercise  the 
filter  and  a  set  of  user-defined  subroutines  which  define  the 
particular  filter  to  be  tested.  A  set  of  input  parameters  is 
varied  for  Input/output  control  and  filter  functioning, 
depending  on  the  data  required  and  the  filter  parameters 
desired  to  be  tested.  SOFE  then  performs  the  number  of  Monte 
Carlo  runs  specified  by  the  user  and  generates  any  output 
files  requested  for  later  evaluation  [61. 

SOFEPL  is  a  plotting  postprocessor  for  SOFE.  when  SOFE 
has  generated  time  history  data  files  for  the  number  of 


specified  Monte  Carlo  runs,  SOFEPL  generates  plot  files  for 
the  types  of  plots  requested,  based  upon  ensemble  averages 
across  all  the  runs  [7]. 


1-^5-.  Organization 

A  more  in-depth  coverage  of  AFTI/F-16  AMAS  concepts  and 
a  brief  presentation  of  the  extended  Kalman  filtering 
techniques  used  in  this  thesis  are  provided  in  Chapter  II. 
Chapter  III  develops  a  constant  turn  rate  (CTR)  EKF  as  a 
baseline  model  for  the  TSE.  A  Gauss-Markov  (GM)  acceleration 
EKF  is  developed  in  Chapter  IV,  as  is  the  U-D  factorization 
implementation  of  the  GM  filter.  Results  of  Monte  Carlo 
simulations  and  an  analysis  of  the  filters  is  accomplished  in 
Chapter  v.  Chapter  V  also  presents  the  conclusions  of  thi3 
study  and  recommendations. 


U_i_  Review  af.  Current 


lu±  Introduction 

This  chapter  presents  a  discussion  of  the  portion  of  the 
AFTI/F-16  aircraft  relative  to  this  thesis.  Including  sensors 
and  the  digital  flight  control  system  (DFCS).  The  concept  of 
integrated  fire/flight  control  (IFFC)  systems  Is  explored, 
and  the  EKF  is  Introduced  as  the  target  state  estimator  (TSE) 
for  the  IFFC.  The  reader  is  expected  to  have  a  background  In 
stochastic  processes  and  Kalman  filters  (estimation  theory); 
for  a  rigorous  development,  see  Reference  4.  This  chapter 
concludes  with  a  discussion  of  the  algorithms  used  in  this 
thesis  to  implement  the  Kalman  filters. 

2^1  AFTI/F-16  Aircraft 

The  AFTI/F-16  is  a  modified  F-16  fighter  aircraft  which 
tests  the  integration  of  new  technologies  into  fighter 
aircraft.  These  modifications  allow  the  AFTI/F-16  to  have 
much  greater  flexibility  than  conventional  aircraft  to 
perform  its  role  as  a  testbed  for  new  technology  (21. 

One  obvious  modification  is  the  addition  of  two  vertical 
canards  to  the  front  underside  of  the  airframe  (see  Figure 
l).  These  vertical  canards  provide  aerodynamic  control 
surfaces  that  allow  the  AFTI/F-16  to  perform  such  maneuvers 
as  direct  side  force  and  vertical  lift  without  pitch 
maneuvering  ( 10 1 . 
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Another  feature  of  the  aircraft  is  its  total  fly-by-wire 
control  system,  which  uses  tr iply-redundant  digital  computers 
to  perform  all  flight  control  and  fire  control  system 
functions.  Changing  system  responses,  aircraft 
configuration,  or  mission  capabilities  is  accomplished  by 
modifying  the  programs  which  run  on  the  system  computers  (2). 


2.2.1  Sensors 

Sensors  for  tracking  a  target  include  the  APG-66  radar 
system,  the  Sensor/Tracker  set  (STS),  and  the  pilot.  For 
purposes  of  this  thesis,  the  pilot  is  not  included  as  a 
sensor  . 

The  APG-66  radar  system  uses  an  integral  Kalman  filter 
to  predict  future  target  position  and  rates.  This  Kalman 
filter  uses  an  inertial  Cartesian  coordinate  frame  to  express 
the  states  of  relative  target  position,  relative  target 
velocity,  and  total  target  acceleration;  i.e.,  the  radar 
states  are 


(2-1) 


where 

B.  is  the  (relative)  target  position  vector, 

¥  is  the  (relative)  target  velocity  vector, 

A  is  the  target  acceleration  vector,  and 
T  indicates  total  (inertial)  quantities 
as  expressed  in  inertial  (e.g.,  North/East/Down)  coordinates. 
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Unfortunately,  target  state  estimates  from  the  radar  are 
considered  to  be  too  noisy  and  Inaccurate  to  be  used  for 
accurate  gun  pointing  [101.  In  fact,  the  least  significant 
bit  of  the  radar  target  position  estimate  represents  16  feet 
(10]  (Note  the  performance  of  the  STS  later  in  this  section). 

The  STS  consists  of  sensors  and  a  digital  computer 
dedicated  for  its  use  to  compute  target  position,  velocity, 
and  acceleration;  i.e.,  as  a  target  state  estimator.  The 
sensors  are  a  FLIR  and  laser  ranger  which  share  common  optics 
to  eliminate  sighting  errors.  The  TSE  currently  employed  in 
the  STS  is  implemented  in  the  same  inertial  Cartesian 
coordinate  frame  as  the  TSE  for  the  radar  and  uses  the  same 
states  as  the  radar  (given  in  Equation  (2-1))  (10). 

Besides  being  used  by  other  on-board  computers  (as 
described  below),  TSE  target  states  are  used  to  generate  rate 
aiding  commands  to  keep  the  FLIR/tracker  head  pointed  at  the 
target  during  target  maneuvers.  Target  position  states  in 
the  TSE  have  an  accuracy  specification  of  two  feet  (one- 
sigma).  This  accuracy  ensures  that  the  fire  control 
computers  can  reliably  predict  future  target  position  [10]. 

When  a  target  is  designated  by  the  pilot,  the  radar 
system's  target  states  and  covariances  are  provided  to  the 
STS,  which  then  begins  tracking.  Data  from  the  INU  are 
available  directly  to  the  TSE,  and  the  TSE  bases  its  updates 
on  a  "snapshot"  of  all  sensor  data.  The  update,  then, 
carries  a  time  tag  so  any  computer  on  the  aircraft  can  tell 


the  precise  time  of  the  most  recent  update.  This  is 


important  because  the  TSE  updates  at  a  30  Hertz  rate  (limited 

,v>. 

by  the  FLIR  scan  rate)  while  all  other  computer  systems  on 
the  aircraft  operate  at  a  50  Hertz  update  rate  (10). 

Digital  Flight  Control  System 
The  digital  computers  on  the  AFTI/F-16  perform  all 
flight  control  and  fire  control  functions.  Figure  II  depicts 
the  interaction  between  various  components  of  the  digital 
flight  control  system  (DFCS).  All  sensor  data  are  digitized 
and  sent  to  one  or  more  digital  computers  for  processing. 

The  flight  control  computers  convert  flight  control  commands 
to  control  surface  commands  so  that  the  aircraft  performs  as 
desired  (101. 


2-JL  Integrated  Elre/-£ll.gnt  control  concepts 

when  weapons  were  first  carried  on  board  aircraft,  the 
pilot  (or  gunner)  was  the  sensor,  data  processor,  and 
controller  for  pointing  and  firing  the  weapon.  Lead  angles 
for  pointing  the  aircraft  or  gun  were  computed  by  the  pilot 
(or  gunner)  during  the  engagement,  and  were  limited  by  his 
experience  and  reaction  time.  To  aid  the  pilot  in  deciding 
when  to  fire  and  how  to  fly  the  aircraft,  target  and  bullet 
prediction  algorithms  such  as  the  lead-computing  sight  and 
the  tracer  algorithm  were  developed.  The  IFFC  aids  the  pilot 
even  more,  by  automating  aircraft  maneuvers  and  gun  control 
commands.  The  aircraft  positions  itself  (via  the  IFFC 
system)  and  fires  the  gun  at  the  target  designated  by  the 
pilot  when  an  acceptable  probability  of  kill  exists  (91. 


Figure  2.  Digital  Flight  Control  System  [10) 


In  order  to  point  the  aircraft  properly,  the  IFFC  system 
computes  bullet  and  target  position  one  bullet  time-of-f 1 ight 
(TOF)  into  the  future.  Corrections  to  the  aircraft  flight 
path  are  then  computed  to  put  the  bullets  on  the  target  at 
the  end  of  the  bullet  TOF.  Future  bullet  positions  can  be 
found  by  using  bullet  trajectory  algorithms  or  by  integrating 
the  bullet  dynamics  equations.  Future  target  position  can 
also  be  determined  by  integrating  the  target  dynamics 
equations  forward  to  the  end  of  the  bullet  TOF  CIO). 

In  order  to  predict  the  future  position  of  the  target 
accurately,  the  IFFC  system  needs  an  accurate  estimate  of 
present  target  states.  In  the  AFTI/F-16,  this  target  state 
estimation  is  accomplished  with  an  EKF  in  the  STS  1101.  This 
EKF  is  the  subject  of  this  thesis. 

ZjJL  Kalman  Filtering 

This  subsection  presents  the  EKF  equations  applicable  to 
this  thesis,  and  illustrates  why  the  conventional  covariance 
updates  are  Inadequate  for  small  wordlength  computers.  An 
alternate  covariance  update  algorithm,  U-D  factorization,  is 
presented  and  used  to  Implement  a  proposed  EKF  for  this 
research.  This  subsection  concludes  with  a  discussion  of  the 
difference  between  the  measurement  frame  and  the  state 


estimate  frame  used  in  the  EKFs  of  this  thesis. 


2.4.1.  Extended  Kalman  Filter  Equations 
The  EKF  equations  for  propagation  and  update  [5:43,44) 
are  presented  below.  For  a  full  derivation  of  these 
equations,  see  references  4  and  5. 

Assume  that  the  system  of  interest  is  described  by  a 
stochastic  process  whose  dynamics  model  is 


£(t) 


f.r£(t),ai.t)/t]  +  G(t)w(t) 


(2- 


where 

£[•,*,•]  is  a  (possibly  nonlinear)  function  of  the 
state  estimates,  inputs,  and  time, 

1£(  * )  1 3  the  estimated  state  n-vector,  and  to)  is 
modeled  as  a  Gaussian  random  n-vector  with  mean  1£o  and 
covariance  Po, 

li(  '  )  is  an  r-vector  of  known  input  functions, 

G ( • )  is  an  n-by-s  measurement  noise  input  matrix, 

and 

!£(•,•)  is  a  zero-mean  white  Gaussian  s-vector 
process  independent  of  £<to)  and  of  strength  Q(t),i.e. 
E(«(  t  )hT(  t+-T) }  =  Q(t)4(t) 

and  the  sampled-data  (discrete-time)  measurement  is  modeled 
as  the  m-vector  process 


i(ti)  =  ll[^(t1),t  1  +  y.(t  )  (2- 

where 

tL(  • ,  *  )  is  a  (possibly  nonlinear)  m-vector  function 
of  the  states  and  time,  and 

y.(  * ,  •  )  is  a  zero-mean  white  Gaussian  m-vector 
process,  independent  of'£(to)  and  *(•,•),  and  of 
covariance  R( t i ) . 

For  propagation  between  measurements,  the  following  two 
equations  are  Integrated  forward  in  time: 


15 


where 


1 1 1  ^  )  ,u( t )  ,  t  ]  (2-4) 

P(tlt.)  =  P  [  t ;“?( 1 1 1  i  )  ]  P  ( 1 1 1  i  )  +  P(tlt.)FT[t/£(tlt.) ] 

+  G(t)Q(t)GT(t)  (2-5) 


F  (  • ,  •  )  i s  defined  to  be  the  n-by-n  matrix  of  partial 
derivatives. 


a  x 


(2-6) 


*=x;(  1 1 1  ) 


For  measurement  update,  the  measurements  z.(t^)  are 
incorporated  by  the  equations 

K  =  P-HT[HP~HT  +  R]-1  (2-7) 

-  tr2:(t-),ti i }  (2-8) 

Pt  =  P-  -  KHP-  (2-9) 

where 


H  is  defined  to  be  the  m-by-n  matrix  of  partial 
derivatives , 


H 


d  hi  2S./ 1  i  J 
d  2L 


t- ) 


(2-10) 


2L(tt)  is  the  measurement  vector  from  sensors  in  actual 
applications  or  a  truth  model  of  the  system  for  simulation 
evaluation. 


and  the  time  arguments  have  been  omitted  for  clarity 
(p-=  P(t^),  etc. ) . 


Equation  (2-7)  calculates  what  is  commonly  called  the 


Kalman  gain  o£  the  filter  and  the  term  in  braces  ({•))  in 
Equation  (2-8)  is  termed  the  residual.  Equation  (2-9) 
performs  the  conventional  covariance  update. 


Note  that  a  matrix  inversion  is  required  to  calculate 

the  Kalman  gain  in  Equation  (2-7).  For  a  matrix  to  be 

invertible,  its  determinant  must  not  be  zero  (the  matrix  must 

not  be  singular).  If  computers  had  infinite  wordlength,  this 

would  not  be  a  problem,  because  the  first  term  (HP“HT)  is 

guaranteed  to  be  at  least  positive  semidefinite  while  the 

second  term  (R)  is  guaranteed  (by  definition)  to  be  positive 

definite  [4:216].  However,  when  a  computer  with  a  finite 

wordlength  (as  all  digital  computers  are)  is  used,  the  filter 

may  diverge  or  totally  fail.  Maybeck  states: 

For  instance,  although  it  is  theoretically  impossible 
for  the  covariance  matrix  to  have  negative  eigenvalues, 
such  a  condition  can,  and  often  does,  result  due  to 
numerical  computation  using  finite  wordlength, 
especially  when  (1)  the  measurements  are  very  accurate 
[eigenvalues  of  R(t^)  are  small  relative  to  those  of 
P(t£),  this  being  accentuated  by  large  eigenvalues  in 
Pol  or  (2)  a  linear  combination  of  state  vector 
components  is  known  with  great  precision  while  other 
combinations  are  nearly  unobservable...  [4:3681. 

In  general,  updating  the  covariance  matrix  usually 

requires  at  least  double  precision  arithmetic  [4:3681.  Even 

double  precision  arithmetic  does  not  guarantee  stability  of 

the  Kalman  recursion  equations;  in  fact,  the  conventional 

Kalman  filter  equations  have  been  shown  to  be  inherently 

numerically  unstable  [1:3381. 

Various  alternate  recursion  relations  have  been 

developed  which  are  inherently  stable.  Notable  among  these 

are  various  square  root  forms  and  U-D  factorization 

algorithms  [41.  The  concept  of  the  square  root  forms  is  to 

propagate  and  update  the  square  root  of  the  covariance 


matrix. 


Such  a  technique  can  yield  twice  the  effective 


precision  of  the  conventional  equations  in  ill-conditioned 
problems,  or  the  same  accuracy  with  about  half  of  the 
wordlength  required  for  the  conventional  equations  [4:3691. 
More  importantly,  the  square  root  forms  have  been  found  to  be 
numerically  stable  [1:3383.  The  U-D  factorization  form 
shares  these  desirable  attributes,  and  has  the  added 
advantage  that  computationally  time-consuming  square  roots 
are  never  explicitly  evaluated  [4:3921. 

2.4.3.  U-D  Factorization 

The  U-D  factorization  method  expresses  the  covariance 
before  and  after  update  at  time  ti  as 

P(t-)  =  U(t“)  D(t-)  UT(t-)  (2-11) 

P(t+)  =  U(t+)  D(t+)  UT(t|)  (2-12) 

where 

U(’)  is  an  upper  diagonal,  unitary  matrix;  i.e.,  all 
elements  below  the  main  diagonal  are  zero  and  all 
elements  on  the  main  diagonal  are  one,  and 

D ( • )  is  a  diagonal  matrix;  i.e.,  the  only  nonzero 
elements  are  on  the  diagonal. 

The  computations  in  Equations  (2-11)  and  (2-12)  are  not 
actually  accomplished;  rather,  algorithms  are  used  to  compute 
s',  U,  and  D  directly  instead  of  and  P. 

Although  the  U  and  D  matrices  are  not  unique,  a  uniquely 
defined  pair  can  be  generated  algor ithmically .  Several 
methods  exist  to  generate  the  n-by-n  U  and  D  matrices;  one 
such  method  (4:3921  is  presented  below: 


First,  for  the  n-th  column 


i  =  n 


!P.  /D 
In  nn 


i  ^  n 


i  =  n 


i  =  n-1,  n-2. 


(2-13) 


Then,  for  the  remaining  columns,  for  j=n-l,  n-2,  ...,  1 


PJJ  :t.  DkkU: 


k  =  jf  1 


l  =  n 


i  J*  n 


i  >  j 


(2-14) 


i  =  j 

P»  1  *  ’-1' 


. 1 


Equations  (2-13)  and  (2-14)  are  normally  used  only  for 
Initialization  of  Uo  and  Do;  propagation  and  update 
algorithms  (such  as  those  below)  compute  U  and  D  factors. 

Measurement  update  using  the  U-D  factorization  filter  is 
first  presented  for  the  scalar  case  [4:394],  then  generalized 
to  the  vector  case  below. 

For  a  scalar  update,  compute 


T  T 

i  =  U  (t-)Hl(ti) 

VJ  =  IVtI,£i 

ao  =  R 


D  1 ,  2 ,  ...,  n 


(2-15) 


Then,  for  k=l,  2,  ....  n,  calculate 


(2-16) 


Dkk(tP  =  °kk(tI)ak-l/ak 
bk  «-  Vk 
pk  =  "fk/ak-l 

,vti)  -  ujk(tr> *  Vk 

b3  «-  bj  *  U3k<CI)vk 


j=l/  2 ,  • • • r  k  1 


where  4-  denotes  replacement,  or  'writing  over'  old 
variables  for  efficiency. 

Vector  measurement  updates  are  processed  component  by 
component,  assuming  the  R  matrix  Is  diagonal  (If  the  R  matrix 
is  not  diagonal,  a  transformation  of  variables  is  required  so 
that  it  is  made  diagonal;  the  diagonal  requirement  is 
discussed  further  below)  (4:3941.  When  U(t+)  and  D(t+)  have 
been  calculated,  the  Kalman  gain  can  be  calculated  from 


K(t . )  =  b/a 
—  l  —  n 


(2-17) 


and  the  states  updated  using  Equation  (2-8)  (4:3941. 

Time  propagation  of  the  U-D  factors  involves  many 
concepts  that  are  beyond  the  scope  of  this  thesis.  One 
method  for  time  propagation  of  U  and  D  is  (4:396-3971: 


Form  the  matrices  Y(t”+^)  and  £(t“+^)  by  augmenting: 


Y(tUl)  =  f^ti+1,tl)U(tt)  I  Gd(tl)l 


(2-18) 


= 


D(tp|  0 
0  I  Q  At.  ) 


(2-19  ) 


where 


l 


** t  i +  1' **  i  *  thS  state  transit^on  matrix 


which  propagates  the  states  forward  in  time  from  t 

t0  ‘i+i' 


is  the  n-by-s  discrete  time  noise  input 


matrix,  and 


Q^j  is  the  s-by-s  noise  matrix  which  has 
statistics  identical  to  Q  at  the  discrete  sample 


times  [4:171] 


vv 


■1. 


1*1 


$(t1+1,T)G(T)Q(T)G  ( 


■)$T  (ti  +  i'T)dT 


(2-20) 


i 


Then,  initialize  n  vectors,  each  of  dimension  (n  +  s), 
through 


[  *1  1  *2 


»n  1  =  Y  (tm> 


(2-21) 


and  iterate  on  the  following  relations  for 
k=n,  n-1,  .  .  . ,  1 : 


^cj^j  ~  C.  j  j  (t  ^  ^  ^  )  a^  ^  ,  j*l,  2,  •••,  n+s) 


Dkk(ti+1)  =  ak£k 


** =  vDkk(trti> 


Ujk(ti+l)  =  S^k 


(2-22) 
3*1,  2,  ...,  k-1 


*3  *3  ’  U3k(fcI+l)ak 


Again,  «-  denotes  'writing  over'  old  variables.  The 
state  estimate  is  given  by 


,  v 


®t-  >  -  *  B . ( t .  )ii( t . ) 


(2-23) 


where 


Bd  is  the  discrete-time  input  matrix,  defined 
by  (4:171] 


1  fcUl 
*'*!• 

«  i. 


/  t )  B  (  t  )  d*r 


( 2-24  ) 


2.1.4. 


Frames : 


Zij.  State 


Whenever  two  quantities  are  summed  they  must  be 
expressed  in  the  same  reference  frame  for  the  sum  to  be 
valid.  However,  in  Equation  (2-7),  the  quantities  utilized 
by  the  STS  are  expressed  in  two  different  coordinate  frames. 
The  TSE  states  and  covariances  are  expressed  in  the  inertial 
reference  frame,  while  the  measurements  and  their  variances 
are  known  only  in  the  LOS  frame.  If  the  frames  are  closely 
aligned  at  a  given  time  (l.e.,  the  sensor  head  is  pointed 


along  the  "north"-axis  of  the  Inertial  frame),  the  variances 
along  each  axis  are  nearly  the  same  in  either  frame  (8:401. 

However,  since  the  reference  frames  are  not  generally 
aligned,  the  covariances  expressed  in  the  LOS  frame  is  skewed 
in  the  inertial  frame  (see  Figure  2.3  for  a  simplified 
example).  Therefore,  a  change  of  reference  frame  must  be 
performed  for  each  update  cycle  (9).  This  change  of  frame 
may  take  one  of  two  forms: 

(1)  The  covariances  and  state  estimates  are  rotated  from 
the  inertial  frame  into  the  measurement  frame  after 
propagation  to  t£  but  before  measurement  update  at  ti>  This 
rotation,  as  used  in  this  thesis,  is  based  upon  an  Euler 


Figure  3.  Covariances  in  Two  Reference  Frames 


angle  transformation  developed  in  Appendix  A  [8].  The 
covariances  and  state  estimates  must  then  be  rotated  back 
into  the  inertial  frame  after  measurement  update  for  normal 
propagation,  a3  in  Equation  (2-5). 

(2)  The  measurement  noise  matrix  R  is  rotated  from  the 
LOS  frame  into  the  Inertial  frame  before  measurement  update 
occurs . 

The  first  option  above  involves  the  use  of  two  matrix 
transformations.  Because  the  transformations  are  Inverse 
(and  transpose)  operations,  calculations  can  be  greatly 
simplified,  as  shown  in  Chapter  III. 

The  second  option  Involves  only  one  transformation,  at 
first  glance.  However,  after  a  diagonal  R  matrix  is  rotated 
into  the  inertial  frame  it  would  not,  in  general,  be 
diagonal.  Since  vector  updates  in  the  U-D  factorization 
algorithm  require  the  R  matrix  to  be  diagonal,  a  change  of 
variables  would  have  to  be  done  to  accommodate  the  update 
algorithm.  These  additional  rotations  are  exactly  the 


,iTr 


rotations  required  In  the  first  option  above;  therefore,  the 
method  shown  in  the  first  option  is  used  in  this  thesis  to 
rotate  the  covariances  and  states  during  update. 


2  Filter  Implementations 

Several  models  have  been  used  in  the  literature  to 
represent  target  motion  in  an  air-to-air  attack  scenario  (81. 
Several  models  use  various  techniques  to  calculate  the 
acceleration  state  as  a  function  of  time  and  Integrate 
successively  to  obtain  the  velocity  and  position.  The  states 
themselves  may  be  expressed  In  a  variety  of  reference  frames. 

One  such  model  Is  the  Gauss-Markov  acceleration  Inertial 
coordinate  (GM)  filter.  In  the  GM  model,  the  accelerations 
are  calculated  using  a  first-order  Gauss-Markov  acceleration 
model,  but  do  not  account  for  persistent  (non-zero  mean) 
accelerations.  Although  the  GM  model  does  not  reflect  the 
real-world  system  with  great  accuracy  (8:17),  the  GM  filter 
has  the  advantage  of  linear  dynamics. 

Another  candidate  model  is  the  constant  inertial  target 
turn  rate  constant-speed  Inertial  coordinate  (CTR)  filter. 

In  the  CTR  filter,  the  accelerations  are  calculated  based 
upon  a  constant  speed  target  going  through  a  constant  turn 
rate,  planar  turn.  The  acceleration  model  of  the  CTR  filter 
approximates  actual  target  accelerations  more  closely  than 
the  GM  filter  [8:3,22],  but  the  dynamics  model  is  nonlinear, 
as  is  shown  in  Chapter  III. 


Because  the  CTR  filter  more  accurately  models  an  aerial 
engagement,  it  Is  used  as  a  baseline  model  for  this  thesis. 
However,  due  to  its  complexity,  the  CTR  filter  is  Impractical 
for  use  in  the  TSE.  The  equations  involved  in  propagating 
and  updating  the  states  and  covariances  require  much  more 
time  and  memory  in  the  CTR  filter  than  in  the  GM  filter,  as 
is  evident  from  the  dynamics  models  of  the  filters  (see 
Chapters  III  and  IV).  For  this  reason,  a  modified  form  of 
the  GM  filter  is  developed  in  Chapter  IV  of  this  thesis  as 
the  proposed  model  for  the  TSE. 


2.5.  Summary 

This  chapter  presented  a  brief  overview  of  AFTI/F-16 
systems  and  IFFC  concepts  to  motivate  target  state  estimator 
Implementation  on  this  aircraft.  Also,  the  Kalman  filtering 
techniques  presented  are  used  as  tools  for  Implementing 


target  state  estimators  in  the  next  two  chapters. 


m 


III ■  constant  Turn  Rate  Filter 

3.1.  Introduction 

The  constant  turn  rate  ( CTR )  EKF  is  presented  in  four 
parts  In  this  chapter.  These  four  parts  are  the  dynamics 
model,  the  measurement  model,  the  noise  models,  and  specific 
implementation  on  SOFE . 

The  CTR  filter  is  a  complex  filter  (as  is  shown  below) 
which  closely  models  the  real-world  system  [8).  Because  the 
CTR  filter  so  closely  models  the  real-world  system.  It  is 
used  as  a  baseline  filter  for  comparison  purposes  and  few 
simplifying  assumptions  are  made  to  reduce  lt3  complexity. 

3.2.  Filter  Dynamics  Modal 

Implementation  of  the  CTR  filter  dynamics  model  Involves 
the  evaluation  of  the  L  and  *  vectors  and  the  G  matrix  of 
Equation  (2-2).  The  jt  vector  (the  estimated  target  state)  is 
defined  by  Equation  (2-1),  expanded  for  clarity: 


r?  1 
1 


/\ 
X  , 


X 


(  3-1  ) 
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where 


*2'  an<^  x3  are  the  relative  target 
positions  along  the  1,  2,  and  3  axes,  respectively,  of 
the  inertial  coordinate  frame  (e.g..  North,  East,  and 
Down ) , 

x4'  ^5'  anc*  ^6  are  the  relative  target 
velocities  along  the  1,  2,  and  3  axes,  respectively,  of 
the  inertial  coordinate  frame,  and 

x?,  $  and  Xg  are  the  total  target 
accelerations  along  the  1,  2,  and  3  axes, 
respectively,  of  the  inertial  coordinate  frame. 

The  differential  equations  describing  the  position 

states  as  a  function  of  the  velocity  states  are  written  as 


(3-2) 

(3-3) 

(3-4) 


The  differential  equations  describing  the  velocity 
states  as  a  function  of  accelerations  are  written  as 


A. 


a  2 
a  3 


(3-5) 

(3-6) 

(3-7) 


where  ai,  aa,  and  a3  are  the  accelerations  of  the 
attacker  along  the  1,  2,  and  3  axes,  respectively,  of 
the  inertial  coordinate  frame. 


The  differential  equations 
states  are  written  as 

*  2  ~ 

x_  =  ( x.  +  v.  ) 

*7  2  -n4  1 

x8  =  -Hwir(x5  +  v2) 

Xg  =  —  Hw.lt  ^  (  Xg  +  Vg) 


describing  the  acceleration 


+  w. 


+  W. 


4-  W. 


(3-8) 

(3-9) 


(  3-10  ) 


where 


Vl'  v2'  anc*  v3  are  the  velocities  of  the  attacker, 
provided  by  the  INO,  along  the  1,  2,  and  3  axes, 
respectively,  of  the  inertial  coordinate  system  (thus,  Xi+v3 
is  the  total  target  velocity  in  the  j-dlrection) , 

W1  •  w2'  an<^  w3  are  zero-mean,  white  Gaussian  noises. 
Independent  of  each  other  and  of  &,  and 

2 

the  term  Hull  is  the  square  of  the  magnitude  of  the 
angular  velocity,  calculated  as 

» 2  .  .2  .  .2 


2  2 
+  A2  *  A3 


(3-11) 


where 
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+ 
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(X, 

+  V 

4- 

(  X, 

)2 

+  ('x, 

(3-12) 
( 3-13) 
( 3-14  ) 
(3-15) 


The  dynamics  equations  for  each  of  the  states  are 
expressed  in  matrix  form  as 


x  - 


S8  -  *2 

x  -  a-> 
9  6 


+  vL) 
-llttH2(x*5  +  v 2  ) 
-ll«ll2(xV  +  v  •  ) 


( 3-16  ) 


which  is  In  the  form  of  Equation  (2-2).  Using  Equation 
(3-16),  the  partial-derivative  matrix  F  is  (see  Equation 
(2-6)  ) 
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F6 

0 

0 

0 

F7 

F8 

F9 

F10 

F11 

F12 

0 

0 

0 

F13 
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F17 
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(  3-17) 


where 


llttH 


2  2($4+Vl)  (A4(A31<,8-A21c‘9)-2(aJ  +k\  )('x,4+v1))' 

+  - ^ - 


(3-18) 


2  ($A  +v,  )  (  AJA,?  -A,xJ  - 


F2  = 


4  1'  4'T9  3  7 


2(A21  .A22  *a\  ><VV  ' 


T 


(3-19) 


2  ( x‘4+v1 )  r  A4(A25?7-A15?8)  -  2  ( A^  +A2  +A3  )(^?6  +  v3)  ) 

T3 - 


F  = 
3 


(  3-20) 


F6  = 


2 ( x^+Vi ) [  A3(x4+Vl)  -  a1(x6+v3)  1 


2(xi+vI)[  A1(x5+v2)  -  A2(x4+v1)  ] 


2(1c,5  +  v2)  (A4(A3'x8-A2x'9)-2(A^  +A2  +A3  )(3?4+v1)] 


( 3-22  ) 


( 3-23) 


(3-24) 


'  .2  2<Vv2>U4<AlVA3V-2(Al  *A2  *A3  "W 

Hull  +  1  —  — —  - 


(3-25) 


F9  = 


2Cx5+v2)[  A4  ( A2'x7-A1Xg  )  -  2  ( A^  +A2  +  Alj  )Cx6+v3)  ] 


2(x5+v2)C  A2{X6+V3}  ~  A3(X5+V2)  1 


3-26) 


(3-27) 


Fll  = 


2  (  X5  +  V2  )  (  A3(x44-v1) 


A1<W  1 


2(‘x5-t-v2)[  A1(‘x5  +  v2)  -  A2(x4+v1)  ] 


(3-28) 


(3-29) 


2(x6  +  v3)(  A4(A3x8-A2x9  )  -  2  ( A^  +a2  fA3  H^+Vj.)  1 


(3-30) 


2 (Xg+v3 ) (  A4 ( A1Xg-A3X7 ) 


2  2  2 
2(A^  +A2  +A3 


(X5+V2)  1 


(3-31) 


•VVjv-V*  %-V  v-' 


)  (x^+V, )  ] 


F18  = 


2(56-v3)(  A 2<Vv3>  -  Wv2>  1 


2($6+v3)(  A3(x4+v1)  -  A1(xg+v3)  ] 


2(x6+v3)(  A1(x‘5  +  v2)  -  A2(x4  +  v1)  ] 


3-32) 


(3-33) 


(3-34) 


(3-35) 


Note  that  the  y.  vector  of  Equation  (2-2),  the  velocity  and 
acceleration  of  the  attacker,  is  imbedded  in  Equations  (3-16) 
and  (3-17).  Also,  the  product  of  G  and  v*  is  given  as  the 
last  (column)  vector  of  Equation  (3-16). 


3.3.  Measurement  Equations 

Four  measurements  are  available  from  the  STS.  Each 
measurement  is  discussed  below  as  it  relates  to  the  states, 
then  the  H  matrix  of  Equation  2-10  is  presented. 


3 , 3 1 Ranae.  Measurement 

A  laser  range  finder  in  the  STS  provides  the  range 
measurement.  The  range,  expressed  in  the  LOS  coordinate 
frame,  is  computed  as  the  time  required  for  a  pul«ed  beam  of 
light  to  travel  to  the  target  and  back,  divided  by  twice  the 
speed  of  light.  The  laser  ranger  in  the  STS  has  a  range 


accuracy  of  2.2  meters  (one-sigma)  (10).  From  spherical 


■ 

i 


3 


geometry,  the  range  to  the  target  is  expressed  as  the  square 
root  of  the  sum  of  the  squares  of  the  relative  position 
states  in  either  the  inertial  or  the  LOS  coordinate  frame; 


i.e.,  using  state  variables, 

R  =  ft*  +  It]  +  V*  )(1/2)  (3-36) 

UJL  Ys.lagl.ty  Measurement 

The  relative  target  velocity  along  the  LOS  is  measured 
by  the  APG-66  radar  system  by  measuring  the  Doppler  shift  of 
the  pulse  of  electromagnet ic  energy  reflected  from  the 
target.  The  accuracy  of  the  velocity  measurement  furnished 
by  the  radar  system  Is  one  foot  per  second  (one-sigma)  CIO). 
The  velocity  measurement  as  a  function  of  the  states  is 
computed,  in  either  the  inertial  or  LOS  coordinate  frame,  by 
forming  the  dot  product  of  the  relative  position  and  relative 
velocity  vectors: 

V  =  (E  •  Y>  /  R  (3-37) 

where 

V  is  the  (signed)  scalar  relative  velocity, 

E.  is  the  relative  position  vector. 


Y  is  the  relative  velocity  vector. 


and  R  is  the  range  to  the  target,  defined  by 
Equation  ( 3-36 ) . 


(3-39) 
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Equation  (3-37)  is  written  using  state  variables  as: 

V  =  (x^  +  5C 5’x'2  +  ”XgX*3 )  /  R  (3-40) 

3 . 3 , J  ,  Angle  Measurements 

The  FLIR  in  the  STS  provides  a  measurement  of  angular 
errors  with  respect  to  the  boresight  of  the  FLIR  head.  The 
boresight,  or  sensor  head  position,  angles  with  respect  to 
the  aircraft  body  are  known  from  position  sensors  in  the  STS, 
and  the  INU  provides  aircraft  attitude  data  directly  to  the 
STS.  Therefore,  all  angles  between  the  inertial  frame  and 
the  LOS  frame  are  known  when  a  sample  occurs,  and  the  azimuth 
and  elevation  angles  from  the  attacker  to  the  target  are 
calculated  as  a  'snapshot'  of  these  angles.  The  STS  measures 
the  azimuth  and  elevation  angles  with  a  precision  of  1.13 
milliradians  (one-sigma)  [101.  This  figure  includes  all 
error  sources:  INU  measurement  errors,  sensor-head  alignment 
errors,  FLIR  pixel  errors,  etc.  [101.  The  angle  measurements 
are  generated  from  the  state  estimates  by  rotating  the  state 
estimates  into  the  sensor  reference  frame  and  comparing 
angular  distances,  as  explained  below. 

Figure  IV  describes  the  angle  relationships  present 
during  a  measurement.  The  Euler  angle  relationships  are 
developed  in  Appendix  A;  the  matrix  which  rotates  quantities 
from  the  inertial  frame  to  the  LOS  frame  is 
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T21  T22  T23 

T31  T32  T33 


(3-41) 


where 


=  cos  0  cos  4  cos  a  cos  0  -  sin  0  cos  4  sin  a 

-  sin  4  cos  8  sin  0  (3-42) 

T^2  =  cos  0  cos  4  sin  b  cos  0  +  sin  0  cos  4  cos  8 

-  sin  4  sin  a  sin  0  (3-43) 

T13  =  -cos  0  cos  4  sin  0  -  sin  4  cos  0  (3-44) 

T21  =  -sin  0  cos  a  cos  0  -  cos  0  sin  a  (3-45) 

T22  =  -sin  ®  sin  *  cos  *  +  cos  •  cos  ®  (3-46) 

T  =  sin  0  sin  a  (3-47) 

23 

=  cos  0  sin  4  cos  a  co3  0  -  sin  0  sin  4  sin  a 

+  cos  4  cos  a  sin  0  (3-48) 


T32  =  cos  0  sin  4  cos  0  +  sin  0  sin  4  cos  a 

+  cos  4  sin  a  sin  0 


(3-49) 


T32  =  -cos  0  sin  4  sin  0  +  cos  4  cos  0  (3-50) 

where 

0  is  the  first  Euler  angle  from  the  inertial  axes 
to  the  body  axes  (the  'azimuth'  of  body  attitude), 

4  is  the  second  Euler  angle  from  the  inertial  axes 
to  the  body  axes  (the  'elevation'  of  body  attitude). 


«  is  the  first  Euler  angle  from  the  body  axes  to 
the  sensor-head  axes  (the  azimuth  angle  of  the  sensor 
head  relative  to  the  body),  and 

0  is  the  second  Euler  angle  from  the  body  axes  to 
the  sensor-head  axes  (the  elevation  angle  of  the  sensor 
head  relative  to  the  body) 

as  defined  in  Appendix  A. 

The  position  state  estimates  are  rotated  into  the  sensor 
frame  by 


and  the  tangents  of  the  azimuth  and  elevation  angles 
calculated  as 


tan 

(a)  = 

& 

X2 

/  -2^ 

(3-52) 

tan 

(e)  = 

& 

*3 

/  ^ 

(3-53) 

Note  that 

there 

is  no 

need  to  correct  the  tangent 

function 

for  a  negative  term  in  the  denominators  of  the  last  two 
equations;  the  denominators  in  either  equation  cannot  be 
negative  if  the  sensor  is  tracking  the  target  within  90 
degrees.  The  tracker  should,  in  fact,  be  pointing  (tracking) 
the  target  within  a  few  tenths  of  degrees,  if  operating 
proper ly . 


3.3.4.  The  Measurement  Matrix 

To  form  the  measurement  vector,  ti(1£/U.,t),  of  Equation 
(2-3),  Equations  (3-36),  (3-40),  (3-52)  and  (3-53)  are 


combined : 


»vy 


Eq . 

(3-36) 

Eq. 

(3-40) 

Eq . 

(3-52) 

Eq . 

(3-53) 

(3-54) 


Evaluation  of  the  partial-derivative  matrix,  H,  of 
Equation  (2-10),  is  straightforward  for  the  range  and 
velocity  measurements.  However,  the  transformation  in  the 
angle  measurements  is  a  function  of  the  position  states,  and 
a  fully  developed  H  matrix  requires  partial  differentiation 
of  terms  of  the  rotation  matrix  T*  with  respect  to  the  state 
vector.  The  terms  involving  differentiation  of  the  rotation 
matrix  are  assumed  to  be  small  compared  to  the  remaining 
terms,  thus  avoiding  this  differentiation  and  the  resulting 
filter  complexity.  This  assumption  has  been  shown  [8:521  to 
be  valid.  The  H  matrix  is  thus  computed  as 


H  = 


H. 


h12  0 


H. 

H. 


H10  HU  0 
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0 

H. 

0 
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0 

Hl 

0 

0 


0 
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0  0 


0 

0 

0 
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( 3-55) 


where 


H,  = 


H, 


H,  = 


/  R 


/  R 


^3  /  R 


=  o2^  -  <35  ^  *5  ^  ^  ]  /  r 


-*\L  x'Jj  t/\L  ✓sLi  ✓sL 


(3-56) 

(3-57) 

(3-58) 

(3-59) 


(3-60) 


H5  =  [ R2x^  -  <2$  +3^  x^  +2^  2!j  )2^  ]  /  R3 

Hg  =  [R2Xg  -  (3^  ^  +3^  x^  )2^  ]  /  R3 


H?  =  2^  /  R 
H8  =  *$  /  R 


Hg  =  X^  /  R 


Hio=  -A  '  (St  >2 


"u*  1  ' 


H12*  ~^3  ^  ^X1  * 


H13*  1  /  “*1 


with 


R  =  [<2*|  )2  +  <x^  )2  +  (2^  ) 2  ]  H/2 ) 

3,1,  Noise  Models 


(3-61) 

(3-62) 

(3-63) 

(3-64) 

(3-65) 

(3-66) 

(3-67) 

(3-68) 

(3-69) 


The  noise  models  used  in  the  CTR  filter  are  divided  into 
two  categories  and  discussed  separately.  The  two  categories 
are  dynamics  driving  noises,  which  comprise  the  &  vector 
(with  strength  matrix  Q)  of  Equation  (2-5),  and  measurement 
noises,  which  comprise  the  y.  vector  (with  covariance  matrix 
R )  of  Equation  ( 2-7 ) . 


3 , 4 JL.  Dynamics  Noises 

The  CTR  filter  acceleration  states  are  driven  by  noises 
w^,  Wj,  and  w^  which  are  assumed  to  be  zero-mean,  white. 


Gaussian,  and  uncorrelated  with  each  other  as  well  as  the 


states.  The  strengths  q^,  q2/  and  q^,  respectively,  of  these 
noises  are  assumed  to  be  equal;  i.e.,  the  target  is  assumed 
to  be  equally  capable  of  acceleration  in  any  direction.  All 
accelerations  are  assumed  to  be  less  than  9  g's,  which  is 
treated  as  a  two-sigma  value  to  account  for  95  percent  of  all 
expected  accelerations.  Thus,  the  strength  of  the  noises 
(r-2)  is  approximately  21,000  feet  2/second  If  the  G  matrix 
of  Equation  (2-2)  is  assumed  to  be  a  9x9  Identity  matrix,  the 
Q  matrix  is 


(3-70) 


where  q1=21000,  1=1,2, 3. 


NQ.1.S&3 


Since  real  sensors  are  imperfect,  the  measurements 
available  from  them  are  corrupted  by  noise.  These  discrete¬ 
time  noises  are  assumed  to  be  zero  mean,  white,  Gaussian,  and 
uncorrelated  with  each  other  as  well  as  the  states.  The  zero 
mean,  white,  and  Gaussian  assumptions  do  an  adequate  job  of 
modeling  modern  sensors,  and  the  noises  are  uncorrelated  if 
the  measurements  are  made  by  Independent  devices.  Although 
the  devices  performing  the  measurements  are  not  totally 


independent,  especially  for  the  angle  measurements,  any 
dependence  is  assumed  to  have  negligible  effects  on  the 
measurements  taken.  If  the  measurements  are  not  independent, 
a  change  of  variables  is  done  (see  Section  2.4.3)  to  effect 
independence.  These  assumptions  yield  a  diagonal  R  matrix, 

’r:  0  0 

0  R  0 

R  ■  0  O  R, 

0  0  0 

where 

Rl'  R2'  R3'  and  R4  are  the  variances  of  the  measurement 
noises  for  the  range,  range  rate,  and  tangents  of  the  azimuth 
and  elevation  angles,  respectively. 

Because  the  standard  deviations  for  the  angle 
measurements  are  small,  the  standard  deviations  for  the 
tangents  of  the  angles  are  approximately  equal  to  the  angular 
standard  deviations  by  the  small  angle  approximation.  The 
measurement  accuracies  presented  in  Section  3.3  are  one-sigma 
values;  the  accuracies  and  variances  are  listed  in  Table  I. 


Table  I 

Measurement  Accuracies  and  variances 


Measurement 

Accuracy  (lx) 

2 

Var lance  (  r  ) 

Range 

6.7  feet 

45  feet2 

Range  rate 

1  foot/sec 

1  feet^'sec2 

Tangent  (az) 

1.13  mrad 

1.2769X10-6rad2 

Tangent  (el) 

1.13  mrad 

1 . 2769X10-6rad  2 

4  J 
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3.5.  C1R  Filter 


Implementing  the  CTR  filter  requires  a  slight 
modification  to  the  general  SOFE  routines  as  well  as 
development  of  filter-specific  subroutines.  A  modification 
to  the  general  SOFE  routines  is  required  so  that  the 
rotations  described  in  Section  2.4.4  can  be  performed  before 
and  after  each  update.  The  filter-specific  subroutines 
define  the  filter  dynamics  model,  the  measurement  model,  and 
the  noise  models,  and  provide  input  of  filter  variables  to 
the  program  (61.  SOFE  contains  many  features  to  enhance  user 
interaction  and  output,  which  are  explained  in  detail  in 
Reference  6. 

SOFE  exercises  the  filter  model  by  integrating  the 
dynamics  equations  forward  in  time  to  the  update  time  via  a 
fifth  order  Kutta-Merson  technique,  then  performing  the 
update  via  a  Carlson  covariance  square  root  algorithm.  The 
integration/update  cycle  is  repeated  until  the  final  time 
specified  by  the  user  is  reached.  The  entire  process  from 
initial  time  to  final  time,  called  a  run,  is  then  repeated  as 
many  times  as  the  user  desires.  Successive  runs  use  a 
different  random  number  sequence  for  the  simulated  input 
noises,  so  averaging  the  outputs  from  a  sufficiently  large 
number  of  runs  can  provide  accurate  performance  sample 
statistics  of  a  certain  filter  configuration  1 6 1 .  This 
ensemble  averaging  is  done  by  SOFEPL  [71.  Results  of  the  CTR 
filter  runs  are  presented  in  Chapter  V.  Before  the  GM  filter 
is  presented,  the  reference  frame  rotation  is  discussed. 


i:fc* 


.» 


2.4.  ,  the  sensor  head  position  with  respect  to  the  inertial 
frame  must  be  known.  In  the  actual  system,  the  3ensor  head 
position  is  given  by  two  Euler  angle  rotations:  The  first, 
from  the  inertial  frame  to  the  aircraft  body  axes,  is 
determined  by  the  INU;  the  second,  from  the  aircraft  body 
axes  to  the  sensor  head  axes,  is  determined  by  angle 
resolvers  in  the  sensor  head  mount  [101.  In  this  thesis  the 
attacker  is  assumed  to  be  benign  (as  discussed  in  Section 
1.3),  so  the  first  Euler  angle  rotation  is  eliminated.  The 
second  set  of  Euler  angles  is  initialized  so  that  the  sensor 
head  is  pointed  directly  at  the  (estimated)  target  position, 
then  incrementally  rotated  at  each  sample  time  by  the 
following  algorithm,  which  is  used  in  the  current 
implementation  of  the  TSE  (101: 

1)  Rotate  the  Inertial  target  position  (from  the 
truth  model,  or  external  trajectory)  into  the  current 


LOS  frame 


P  s  t  P 
-L  *L  -I 


( 3-72  ) 


where 


p  is  the  target  position  vector  in  the  LOS  frame, 
L 

T*  is  the  Euler  rotation  matrix.  Equation  3-41, 


P.  is  the  target  position  vector,  inertial  frame. 


2)  Find  the  error  angles  between  LOS  and  the 


target  position  vector: 


ea  =  arctan  (£2l  /  £1L) 


( 3-73  ) 


efi  =  arctan  [P^  /  (P2L  +  P^L  )(1/2))  (3-74) 

where 

ea  az^mu^h  pointing  error  angle/ 

ee  is  the  elevation  pointing  error  angle,  and 
PiL  is  the  component  of  the  PL  vector  in  the  i- 
direction,  1=1, 2,3. 

3)  Calculate  the  velocity  of  the  target  normal  to 
the  range  vector  by  the  cross  product  operation 

ar  ■  Ej  X  Xj  (3-75) 

where 

B-j  is  the  (estimated)  target  velocity  normal  to 
the  range  vector,  expressed  in  the  inertial  frame, 

is  the  (estimated)  Inertial  target  position, 

X  is  the  (estimated)  inertial  velocity  vector, 
and  A 

X  indicates  the  cross  product  operation. 

4)  Rotate  the  velocity  normal  to  the  range  into 
the  LOS  frame,  and  convert  to  unit  vectors 

*L  =  Tl  Rj  /  MBjM2  (3-76) 

where 

VJL  is  the  rate  aiding  vector  in  the  LOS  frame,  and 
2 

HRjII  is  the  square  of  the  relative  target  range. 

5)  Calculate  the  rate  aiding  terms  as 

W  (1)  =  Wo(l)  +  26 . 46 • At  *  e  (3-77) 

n  a 

W  (2)  =  Wo(  2  )  +  2  6 . 4  6 • At  *  e  (3-78) 

n  e 


where 


Wn ( i )  is  the  new  rate  aiding  term  in  the  i- 
direction,  i=l,2,3, 

Wo(i)  is  the  rate  aiding  term  in  the  i-direction 
from  the  previous  sample  time  (Initialized  to  zero),  and 

At  is  the  sample  period  (time  between  updates). 


6)  Calculate  incremental  gimbal  rotations 
next  propagation  cycle  as 

Aa  =  (^(3)  +  12.6-ea  +  Wn(l))At 
Ae  =  (5^(2)  +  12 . 6  •  ee  +  W  (2))At 

where 

Aa  is  the  amount  the  gimbal  rotates  in  azimuth  in 
the  next  At  time  interval, 

Ae  is  the  amount  the  gimbal  rotates  in  elevation 
in  the  next  At  time  interval,  and 

4^(1)  is  the  i-th  component  of 


over  the 

( 3-79  ) 
(3-80) 


7)  Calculate  the  Euler  rotation  matrix  TAt  which 
rotates  the  sensor  head  to  its  desired  position  at  the 
next  sample  time,  based  upon  the  incremental  gimbal 
rotations  from  Equation  3-41. 


8)  Multiply  the  body-to-LOS  Euler  rotation  matrix 
by  the  incremental  gimbal  rotation  matrix  to  yield  the 
body-to-LOS  Euler  rotation  matrix  for  the  next  update: 


where 


(  3-81) 


next 


,B 


is  the 


Ln 

update  time. 


is  the 


body-to-LOS  rotation  matrix  for  the 
incremental  gimbal  rotation  matrix,  and 


A  f- 


TLq  is  the  body-to-LOS  rotation  matrix  for  the 
present  (just  completed)  update. 

The  transformation  from  inertial  to  LOS  is  given  at  the 
next  update  time  by  multiplying  the  inert ia 1-to-body  rotation 
matrix,  obtained  from  attitude  information  from  the  INU,  by 
the  new  body-to-LOS  rotation  matrix.  Since  a  benign  attacker 
is  assumed  in  this  thesis,  the  Inertial-to-body 
transformation  matrix  (for  this  research  only)  is  just  a 
three-by-three  identity  matrix. 


3.6.  Coordinate  Transformations 

As  discussed  in  Section  2.4.4,  the  covariances  and 
states  must  be  rotated  into  the  LOS  frame  for  update  and  back 
into  the  inertial  frame  for  propagation.  The  states  are 
rotated  into  the  LOS  frame  by  the  transformation  matrix  given 


Equation  3-41 

before  update: 

- 

TL  °3X3 

0  3X3 

*L  = 

°3X3  TL 

°3X3 

*1 

(  3-82 

0  3X3  0  3X3 

P 

t-i  j 

E- 

where 

is  the  state  estimate  before  update  at  t^  in  the  LOS 

frame , 

is  the  state  estimate  before  update  at  t^  in  the 
inertial  frame,  and 

t£  is  the  inertial-to-LOS  transformation  matrix 
computed  by  Equation  (3-41). 


Since  the  transformation  matrix  TL  used  for  this  rotation  is 
orthogonal.  Equation  (3-82)  can  be  used  to  transform  the 
states  from  the  LOS  frame  to  the  inertial  frame  after  update, 
if  the  transpose  of  T*  is  used  in  place  of  T*  and  the  roles 
of  and  }?r  are  interchanged.  The  covariances  are  rotated 

I  Li 

into  the  LOS  frame  in  much  the  same  manner,  except  a  post- 
multiply  by  the  transformation  matrix  is  required: 

PL  = 


°3X3 

°3X3 

tl 

0  3X3 

°3X3 

°3X3 

Tl 

°3X3 

p; 

°3X3 

°3X3 

°3X3 

°3X3 

0  3X3 

0  3X3 

7l 

where 

P£  is  the  covariance  matrix  before  update  at  in  the 
LOS  frame, 

P~  is  the  covariance  matrix  before  update  at  in  the 
Inertial  frame,  and 

Tjf  is  the  inert ial -to-LOS  transformation  matrix 
computed  by  Equation  (3-41). 

3.7.  Summary 

The  CTR  filter  is  presented  in  this  section.  The  filter 
dynamics  equations,  measurement  equations,  and  noise  models 
are  presented,  and  implementation  on  SOFE  is  discussed. 
Finally,  the  method  used  in  this  thesis  to  rotate  the  states 
and  covariances  between  the  Inertial  and  LOS  coordinate 
frames  is  presented. 

Because  it  models  the  real-world  situation  well,  the  CTR 
filter  is  used  as  a  baseline  in  this  thesis  and  few 


assumptions  are  made  to  simplify  its  complexity.  The  next 
chapter  presents  a  simpler  filter  model  and  makes  more 
simplifying  assumptions  to  reduce  recursion  computation  time. 
Results  from  both  filter  models  are  compared  in  Chapter  V. 
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12Ll. 


Filter 


4-.JL.  introduction 

Two  implementations  of  the  Gauss  Markov  (GM)  EKF  are 
presented  in  this  section.  The  first  implementation  uses  the 
propagation  and  update  techniques  presented  in  Section  3.5 
for  the  CTR  filter.  The  second  implementation  uses  the  U-D 
covariance  factorization  techniques  described  in  Section 
2.4.3  for  propagation  and  update.  Both  implementations  are 
shown  to  be  less  complex  than  the  CTR  EKF  presented  in 
Chapter  III.  Reduced  complexity  is  desirable  because  the  TSE 
on  board  the  AFTI/F-16  is  limited  in  both  memory  allocation 
(32K  of  program  space  and  32K  of  dynamic  RAM)  and  processing 
time  (7  milliseconds  per  update  frame)  [10] . 

4.2L,  F-ilter.  Dynamics  Model 

The  dynamics  model  of  the  GM  filter  uses  the  same  states 
as  the  CTR  filter,  given  in  Equation  (3-1).  The  dynamics 
model  for  the  GM  filter  is 


where 


£ 


vl'  w2'  an<^  w3  are  zero-mean,  white  Gaussian 
noises  independent  of  each  other  and  of 

Tl'  t2'  an<*  t3  are  correlation  time  constants 
for  the  acceleration  states,  and 

al'  a2'  and  a3  are  the  attacker  accelerations  in 
the  1,  2,  and  3  directions. 

The  correlation  time  constants  characterize  the  half¬ 
power  point  of  the  power  spectral  densities  of  the  correlated 
acceleration.  The  frequency  at  which  the  half-power  point 
occurs  is  assumed  to  be  at  about  two  Hertz  for  high- 
performance  aircraft  [8:613,  yielding  correlation  time 
constants  of  two  seconds.  The  correlation  time  constants  are 
assumed  to  be  equal  in  all  three  directions  because  the 
target  is  assumed  to  have  the  same  acceleration  capacity 
along  any  axis  of  the  Inertial  frame  [8:613. 

Since  the  dynamics  model  is  linear  and  time  invariant, 
the  dynamics  matrix  F  of  Equation  (2-6)  is  calculated  as 
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(  4-2  ) 


which  Is  obviously  a  much  simpler  matrix  to  evaluate  than 
the  F  matrix  for  the  CTR  filter  (Equation  (3-17)). 


Measurement  Equations 

In  the  GM  filter,  the  states  and  covariances  are 
propagated  and  updated  in  the  inertial  frame,  so  the 
measurements  are  accomplished  differently  than  in  the  CTR 
filter.  The  range  and  range  rate  measurement  equations  are 
are  identical  in  any  frame,  and  the  covariances  associated 
these  measurements  are  equal  in  either  the  inertial  or  LOS 
frame.  However,  the  covariances  of  the  angle  measurements  are 
angles  about  the  LOS  vector,  which  are  not  the  same  in  the 
different  frames  unless  they  describe  a  circular  cone  about 
the  LOS  vector.  Unfortunately,  the  covariances  describe  a 
squared  cone,  as  shown  in  Figure  5.  The  covariances  are 
equal  if  the  planes  described  by  the  1-2  axes  of  the 
coordinate  frames  (the  'horizontal'  planes)  are  aligned,  and 
are  skewed  up  to  45  degrees  otherwise.  The  FLIR  head  is 
roll-stabilized  in  the  aircraft  (10),  30  the  'horizontal' 
planes  should  always  be  closely  aligned. 


Figure  5.  Angular  Error  Covariances 


Since  the  states  are  not  rotated  Into  the  LOS  £rame  for 


update,  different  equations  are  used  for  the  angular 
measurements:  the  angles  are  measured  as  Euler  angles;  i.e.. 


tan  (a)  =  1?2  /  (  4-3  ) 

tan  (e)  =  ^  )(1/2)  (4-4) 

and  the  range  and  range  rate  measurements  are  given  in 
Equations  (3-36)  and  (3-40),  respectively. 

The  H  matrix  is  calculated  as 


hl  h2  h3  0  0  0  0 

H4  H5  H6  H1  H2  H3  0 

H7  H0  0  0  0  0  0 

H9  H10  H11  0  0  0  0 
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(4-6) 

(4-7) 

(4-8) 

(4-9) 

(4-10) 

(4-11) 

(4-12) 

(4-13) 

(4-14) 

(4-15) 

(4-16) 


4.4.  Noise  Models 


As  with  the  CTR  filter,  two  types  of  noise  models  are 
apparent  in  the  GM  filter:  measurement  and  dynamics  noises. 
Measurement  noise  models  are  identical  for  either  filter 
because  the  same  devices  are  used  to  perform  the 
measurements.  Dynamics  noise  models  are  acquired  by  assuming 
that  zero-mean,  white,  Gaussian  noises  drive  the  acceleration 
states  and  performing  a  steady-state  covariance  analysis. 
According  to  Worsley  {8:631,  the  results  of  this  covariance 
analysis  results  in  the  equation 

2  P33/t  =  qx  (4-18) 

where 

P33  is  the  covariance  associated  with  the 
acceleration  state  along  any  axis,  and 

q.^  is  the  strength  of  the  driving  noise. 

Assuming  the  same  value  of  error  in  the  target  acceleration 

estimate  as  in  the  CTR  filter  and  the  correlation  time 

constant  of  2  seconds,  the  value  of  q^,  1=1, 2, 3  is  21000 
2  5 

feet  per  second  ,  the  same  as  for  the  CTR  filter. 

4,5-  Acceleration  Rotation 

The  CTR  filter  dynamics  equations  use  a  constant  target 
turn  rate  model,  as  developed  in  Appendix  B,  to  correlate  the 
acceleration  states  with  the  velocity  states  in  an  effort  to 
match  the  expected  target  performance  in  an  aerial  engagement 
[8:3,221.  The  GM  filter  dynamics  model  does  not  correlate 


the  acceleration  states  to  anything,  but  rather  assumes  that 
the  target  is  always  capable  of  equal  accelerations  in  all 
directions.  Thus,  the  GM  filter  does  not  effectively  model  a 
real  aircraft,  because  most  aircraft  accelerations  occur 
normal  to  the  aircraft's  flight  path.  One  method  of  forcing 
the  GM  filter  to  model  a  real  aircraft  more  closely  is  to 
rotate  the  acceleration  states  to  make  them  perpendicular  to 
the  velocity  states  just  before  the  states  are  updated  (10). 
The  acceleration  rotation  is  accomplished  using  two  vector 
cross  products,  as  shown  below: 

^  *  £j,  X  (  A0  X  £j,)  (4-19) 

where 

is  the  new  (rotated)  acceleration  vector, 

AQ  is  the  old  ( non-rotated )  acceleration  vector, 

^  is  the  target  total  velocity  vector,  and 
X  denotes  the  cross  product  operation. 

For  both  implementations  of  the  GM  filter,  the  rotation 
above  is  performed  on  the  state  estimates  after  propagation 
but  before  update.  Thus,  the  GM  filters  used  in  this  thesis 
are  modified  GM  filters. 

1.1 6  .  SM  Filter.  Implementation 

As  mentioned  in  Section  4.1,  the  GM  filter  is 
implemented  two  ways  in  this  thesis.  The  first 
implementation  mirrors  the  CTR  filter  implementation,  shown 
in  Section  3.5,  using  Carlson  square  root  covariance  updates 
and  propagation  by  Integration.  The  second  implementation 


uses  U-D  covariance  factor ization  techniques  to  update  and 
propagate  the  covariance  matrix,  and  the  state  transition 
matrix  to  propagate  the  filter  states.  An  overview  of  these 
techniques  is  presented  in  Section  2.4.3;  the  remainder  of 
this  section  is  devoted  to  explaining  the  implementation  of 
these  techniques  using  SOFE. 

The  state  transition  matrix,  d>(t,to),  propagates  the 
states  from  to  to  t  and  is  used  to  propagate  the  covariances. 
For  a  linear,  t ime- invar iant  system  dynamics  model  such  as 
the  GM  filter,  the  state  transition  matrix  can  be  calculated 
as  an  inverse  Laplace  transform  [4:421: 

<t>(t-to)  =  £*1{[sl  -  Fl_1i  (4-19) 

where 

£  1  is  the  inverse  Laplace  transform  operator,  and 
s  is  the  Laplace  integration  variable. 


The  state  transition  matrix  for  the  GM  filter  is 
calculated  (assuming  r±=r2=r3=r)  as 
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(4-21) 


54 


where 

2  2 

=  T  At  -  T  +  r  exp  ( -At/ T  ) 
$2  =  t  -  ▼  exp(-At/T) 

$3  =  exp(-At/T) 


(4-22) 

(4-23) 

(4-24) 


Since  Equations  (4-22)  and  (4-23)  are  differences  of  numbers 
of  similar  magnitude,  their  calculation  may  be  inaccurate  if 
done  by  a  small  wordlength  computer,  especially  for  small 
values  of  At.  For  this  reason.  Equations  (4-22)  to  (4-24) 
are  expanded  in  a  Taylor's  series.  The  series  is  trucated 
after  the  first  two  terms,  resulting  in  the  approximations 
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3  n=0  n!  Tn 

These  approximations  introduce  a  maximum  error  of  0.0014 
percent  into  the  calculations  for  5. 

To  propagate  the  states  if  an  input  (attacker 
acceleration,  see  Equation  (4-1))  were  present,  the  Bd  matrix 
of  Equation  (2-23),  as  defined  in  Equation  (2-24),  is 
required  (4:171): 
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(  4-28  ) 


where  B  Is  the  continuous-time  input  matrix  imbedded  in 
Equation  (4-1).  Assuming  the  input  is  the  three-element 
attacker  acceleration  vector,  the  d i3crete-t ime  input  matrix 


is  calculated  as 


AtV2  0 
0  At2/2 


0  AtV2 
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(4-29  ) 


where  At  =  t^+1  -  t^ 

The  covariances  are  propagated  by  using  Equations  (2-18) 
to  (2-22).  The  values  for  and  in  these  equations  is 
calculated  by  assuming  that  G,  the  continuous-time  input 
matrix,  is  a  nine-by-nine  identity  matrix  and  Q,  the 
continuous-time  noise  matrix,  is  given  by  Equation  (3-70). 
Then,  Equation  (2-20)  is  used  to  calculate  [4:171]  as 
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(4-33) 
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The  errors  incurred  by  truncating  the  series  after  the  first 
two  terms  are  less  than  0.0094  percent  of  each  Q  value. 

In  the  equations  for  propagating  the  covariances 
(Equations  (2-18)  to  (2-22)),  the  matrix  is  assumed 
diagonal  and  the  matrix  is  upper  triangular.  The  required 
forms  of  these  matrices  are  acquired  by  breaking  the 
matrix  of  Equation  (4-30)  into  an  upper  diagonal  unitary 
matrix  and  a  diagonal  matrix  using  the  algorithm  presented  in 
Equations  (2-13)  and  (2-14).  Note  that  this  is  the  same 
algorithm  used  to  generate  the  U  and  D  factors  of  the  initial 
covariance  matrix,  so  additional  programming  is  not  required. 

To  Implement  the  second  form  of  the  GM  filter  in  SOFE, 
several  modifications  are  required.  First,  the  integration 
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routines  are  removed  and  replaced  with  the  propagation 
algorithms  presented  in  Section  2.4.3.  Then,  the  Carlson 
square  root  update  algorithms  are  replaced  with  the  U-D 
update  algorithms,  also  presented  in  Section  2.4.3.  The 
filter  dynamics  models  are  removed,  because  these  models  are 
imbedded  in  the  propagation  algorithms.  Much  of  the 
remainder  of  the  program  is  then  modified  to  allow  several 
input/output  calls,  embedded  in  the  removed  routines,  to  be 
performed . 

4  ,.l.  Summary 

The  GM  EKF  is  presented  in  this  chapter.  The  GM 
dynamics  model  is  shown  to  be  less  complex  than  the  CTR 
filter  dynamics  model  of  the  previous  chapter.  Also,  by  not 
rotating  the  states  and  modifying  the  measurements,  the 
filter  becomes  less  complex  to  simulate  because  the  sensor 
head  need  not  be  manipulated  in  order  to  acquire  the 
inertial-to-LOS  transformation  matrix. 

In  Chapter  V,  results  of  the  GM  EKFs  are  compared  to  the 
results  from  the  CTR  filter  of  Chapter  III. 


5.1. 


This  chapter  presents  the  results  of  the  SOFE 


simulations  of  the  filters  developed  in  Chapters  r I I  and  IV. 
Results  from  each  filter  are  addressed,  then  the  tracking 
capabilities  of  all  the  filters  are  compared.  Finally, 
conclusions  of  the  research  and  recommendations  for  further 
study  are  presented. 

Each  SOFE  simulation  exercised  a  filter  over  20  Monte 
Carlo  runs  to  monitor  the  statistical  performance  of  the 
filter.  Each  of  these  runs  is  made  against  a  file  of 
trajectory  data  that  is  Invariant  from  run  to  run  and  filter 
to  filter.  In  this  way,  the  relative  performance  of  the 
various  filters  can  be  compared. 


5.X-2.  Trajectory  G&naratlcn 

Trajectories  for  the  simulations  are  generated  by 
computing  the  effects  of  an  acceleration  on  a  point  mass 
traveling  at  a  specified  speed.  Any  changes  in  applied 
accelerations  were  accomplished  at  the  rate  of  nine  g's  per 
second.  Only  basic  planar  moves  were  attempted;  the  point 
mass  first  executed  a  nine  g  right  turn  at  three  seconds  into 
the  simulation,  then  a  nine  g  left  turn  at  five  seconds.  At 
eight  seconds,  the  horizontal  accelerations  dropped  (over  a 
one-second  interval)  to  zero.  Then,  a  similar  move  was  made 
in  the  vertical  plane;  first,  a  5  g  downward  turn  at  9 
seconds,  then  a  5  g  upward  turn  at  11  seconds,  and  falling 


y  /  /  . -  .*  /  /  /  y 


/.  ,*■  \  \ , 


off  to  0  at  14  seconds.  The  target  is  thus  propagated 
forward  from  a  given  set  of  Initial  conditions  and  compared 
to  the  position  of  the  attacker,  which  is  always  headed  due 
north  at  1000  feet  per  second.  The  relative  position, 
relative  velocity,  and  total  target  position,  all  referenced 
to  the  Inertial  coordinate  frame,  are  then  stored  in  a 
trjectory  file. 

Two  sets  of  initial  conditions  generated  two  different 
trajectories.  The  first  set  of  initial  conditions  places  the 
target  5000  feet  directly  ahead  of  the  attacker,  going  in  the 
same  direction  at  the  same  speed,  as  depicted  in  Figure  VI. 
The  second  set  of  initial  conditions  places  the  target  10,000 
feet  away  at  a  -45  degree  angle  from  the  heading  of  the 
attacker  (to  the  attacker's  left).  In  the  second  trajectory, 
the  target  is  travelling  at  1100  feet  per  second  and  at  a  45 
degree  heading;  i.e.,  the  target  is  crossing  in  front  of  the 
attacker,  and  the  first  accelerations  are  applied  towards  the 
attacker  (see  Figure  VII). 

The  trajectories  are  computed  at  a  3000  Hertz  rate  to 
ensure  high  accuracy  of  the  simulation  results.  However, 
data  are  output  to  the  trajectory  file  at  a  30  Hertz 
(simulated)  rate  to  be  compatible  with  the  filter  sample 
rate.  If  the  data  rate  of  the  trajectory  file  is  different 
from  the  system  sample  rate,  SOFE  Interpolates  the  trajectory 
data  as  needed  using  cubic  splines.  This  interpolation 
requires  significant  processing  time  each  cycle  (frame)  [61. 


LJL-  Filter 


Each  filter  is  run  on  SOFE  for  20  Monte  Carlo  runs  for 
each  trajectory  file.  Twenty  Monte  Carlo  runs  is  considered 
to  be  adequate  based  upon  an  analysis  done  by  Worsley  [81. 

For  each  filter/trajectory  combination,  an  ensemble  average 
across  all  20  runs  is  obtained.  Plots  of  these  ensemble 
averages  are  presented  in  Appendix  C.  The  plots  reflect  the 
error  between  the  filter  estimated  state  and  the  trajectory 
state,  and  the  covariance  associated  with  each  state,  one 
state  per  plot.  An  analysis  of  the  plots  follows;  Table  II 
summarizes  the  maximum  position  error  along  each  axis  and  the 
maximum  total  position  error  for  each  of  the  filter/ 
trajectory  combinations. 

The  CTR  filter  keeps  the  error  between  the  estimated 
state  and  the  trajectory  data  well  within  the  covariance 
boundaries  for  all  but  a  short  time  during  the  simulation. 
Also,  the  position  states  always  have  less  than  three  feet 
total  error  along  any  axis.  The  CTR  filter  appears  to  be 
tracking  well. 

Two  different  results  are  obtained  for  the  two 
Implementations  of  the  GM  filter.  The  filter  which 
Integrates  the  equations  forward  in  time  had  results  nearly 
Identical  to  those  of  the  CTR  filter.  However,  the  GM  filter 
which  used  propagation  and  U-D  covariance  factorization 
techniques  tracked  comparat ively  poorly;  the  maximum  total 


(UD) 
CTR 
GM 

(UD) 


Apparently,  a  flaw  exists  In  the  Implementation  of  the 
U-D  factorization  GM  filter.  The  U-D  filter  13  tracking,  but 
not  quite  as  well  as  expected,  given  the  performance  of  the 
other  two  filters.  Varying  system  parameters  seems  to  make 
performance  worse  rather  than  better.  Testing  of  the  various 
routines  reveals  that  the  flaw  is  probably  in  the  routines 
which  propagate  the  states  and  covariances  forward  in  time, 
and  not  in  the  update  routines.  A  possibility  exists  that 
SOFE  does  not  work  properly  without  some  of  the  routines 
which  are  removed  for  this  filter  implementation. 

Not  rotating  the  states  and  covariances  into  the  LOS 
frame  for  update  seems  to  have  very  little  effect  on  filter 
performance.  The  CTR  filter,  which  did  the  rotation, 
performs  little  better  than  the  Integrated- forward  GM  filter. 


which  does  not  do  the  rotation. 


LjJLl.  Recommendations 

Each  of  the  filters  implemented  in  this  thesis  uses 
numerically  stable  covariance  update  techniques.  Two  quite 
different  dynamics  models  are  used  with  comparable  results, 
and  the  covariance  matrix  rotation  during  update  has  no 
significant  effect  on  filter  performance.  However,  the 
previous  AFTI/F-16  TSE  does  not  perform  as  desired  [21. 

Since  the  filter  dynamics  models  and  covariance  rotations  had 
no  significant  effect,  the  effect  of  using  conventional 
covariance  update  techniques  should  be  investigated. 

Also,  the  effects  of  finite  wordlength  should  be 
investigated  in  more  detail.  The  equations  developed  in  this 
thesis  can  be  run  on  a  comparatively  short  wordlength 
computer,  such  as  a  Zenith  Z -248,  to  study  the  effects  of 
truncated  wordlength3,  or  the  wordlength  of  each  quantity  can 
be  trucated  after  each  operation  to  simulate  the  shorter 
wordlength . 
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Ll  Euler  Anql.ft 


This  appendix  briefly  presents  the  Euler  angle 
relationships  used  in  this  thesis.  Three  reference  frames 
are  used:  the  inertial  frame,  which  is  stationary  for  all 
time,  the  (attacker  aircraft)  body  axis,  which  may  both 
translate  and  rotate,  and  the  tracker  (or  LOS)  frame. 

The  first  reference  frame,  the  inertial  frame,  is 
defined  by  the  orientation  of  the  INU  and  denoted  by  a 
subscript  'I'.  The  origin  of  the  axes  is  defined  as  the 
center  of  the  INU,  because  all  rotations  are  measured  around 
that  point.  A  common  name  for  this  type  of  inertial  frame  is 
a  wander  azimuth  frame,  because  the  'down'  direction  is 
always  defined  and  'ahead'  or  'right'  can  be  defined  in  any 
convenient  manner. 

The  second  reference  frame,  the  body  frame,  is  defined 
to  have  the  1-axis  pointed  out  the  nose  of  the  aircraft,  the 

2- axis  pointed  along  the  right  wing,  and  the  3-axis  pointed 
out  the  'bottom'  of  the  aircraft  body.  The  body  frame, 
denoted  by  the  subscript  'B',  has  its  center  at  the  center  of 
gravity  of  the  aircraft. 

The  third  reference  frame,  the  LOS  frame,  is  defined  by 
the  orientation  of  the  sensor  head.  The  1-axis  points 
perpendicular  to  the  plane  of  the  sensor  (FLIR)  array  and  the 

3- axis  points  'down'  (the  sensor  head  is  r ol 1 -s tabi  1 ized ) . 

The  2-axis  completes  the  right-handed  coordinate  frame  (it 
points  to  the  'right'). 
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Each  reference  frame  can  be  related  to  the  next  by  Euler 


angle  rotations.  These  rotations  are  defined  by  two  angles, 
say  9  and  These  angles  define  how  much  the  base  reference 

frame  is  rotated,  first  to  the  right  in  azimuth  then  down  in 
elevation,  to  match  the  orientation  of  the  second  reference 
frame . 

Although  Euler  angles  are  easily  defined  for  most 
relative  orientations,  they  are  ill-conditioned  if  the  second 
Euler  angle  approaches  90  degrees.  If  the  second  Euler  angle 
equals  90  degrees,  the  first  Euler  angle  is  undefined  because 
the  1-axis  of  the  second  reference  frame  cannot  be  projected 
into  the  1-2  plane  of  the  base  reference  frame.  The 
possibility  of  this  occurring  is  assumed  to  be  negligible  in 
this  thesis;  other  coordinate  trans f ormat 1 ons ,  such  as 
direction  cosines  or  quaternions,  could  be  used  to  eliminate 
any  problem  with  reference  frame  orientations. 

However,  because  Euler  rotations  are  performed, 
especially  in  the  INU  where  the  STS  has  no  control,  one 
property  of  Euler  rotations  is  especially  important  in  this 
thesis.  This  property  is  the  following:  if  two  or  more  Euler 
rotations  are  performed  to  go  from  the  basis  reference  frame 
to  the  measurement  reference  frame,  the  orientation  of  the 
measurement  reference  frame  with  respect  to  the  basis 
reference  frame  can  not .  in  general,  be  described  by  a  single 
Euler  rotation  using  the  3ums  of  the  rotation  angles. 
Therefore,  because  the  angle  measurements  are  performed  in 
the  sensor  head  frame,  the  position  of  the  sensor  head  must 
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be  known  (or  estimated)  in  the  Kalman  filter.  Estimated 
angles  are  then  rotated  into  the  sensor  head  reference  frame 
and  compared  to  the  actual  measurement.  Conversely,  the 
measurement  could  be  rotated  into  the  basis  (inertial)  frame 
and  compared  to  the  estimated  line-of-sight;  this  thesis  uses 
the  former  approach  because  the  covariance  matrix  undergoes 
the  same  rotations  during  measurement  update  (as  explained  in 
Chapter  II). 


2 

The  coefficient  Hull  used  in  the  constant  turn  rate 

( CTR )  filter  is  derived  in  Reference  8  and  included  here  to 

make  this  thesis  as  complete  a  work  as  possible. 

2 

The  CTR  coefficient  Hull  ,  which  is  the  square  of  the 
magnitude  of  the  target's  inertial  turn  rate,  is  developed 
from  the  application  of  the  Coriolis  theorem,  written  as 

Zp  +  tt1  X  Zp  (B-l) 

where 

Zp  1 s  the  inertial  target  velocity, 

it1  is  the  inertial  target  angular  velocity, 

X  denotes  the  cross-product  operation,  and 

the  superscripts  I  and  T  before  the  derivatives  indicate 
that  the  derivatives  are  taken  in  the  inertial  and  target 
body  frame,  respectively. 

Now,  the  first  term  of  the  right-hand  side  of  Equation  (B-l) 
is  zero  because  the  target  is  assumed  to  be  at  constant 
speed.  Taking  the  derivative  of  the  remainder  of  Equation 
(B-l)  with  respect  to  time  yields 

1  j  2  I  j 

d  d 

-  Zp  =  —  (fcl  X  5^  )  (B-2) 

dt2  dt 

or,  expressed  in  the  target's  body  frame. 


Now,  since  both  the  target  speed  and  angular  velocity  are 
assumed  constant,  the  first  term  on  the  right-hand  side  of 
Equation  ( B - 3 )  is  zero,  yielding 


1,2 

^  II 

-  Zj,  =  * X  (* 1  X  Jtp)  (B- 

dt2 

Using  the  relationship  for  a  triple  cross  product.  Equation 
(B-4)  is  written  as 

1,2 

^  I  III 

-  Zj,  =  •  Z£)  1*  ~  •  *  )  Zj,  (B- 

dt2 

The  first  term  of  Equation  ( B  —  5 )  is  zero  since,  for  a 
planar  turn,  the  target's  inertial  velocity  and  angular 
velocity  vectors  are  perpendicular.  Thus,  Equation  (B-5) 
becomes 

ld2 

-  Zv  =  -  Httll2  ZT  (B- 

dt2 

which  is  the  vector  form  of  the  derivatives  of  the 

acceleration  states  of  the  CTR  filter. 

2 

To  compute  If mH  ,  tne  target's  inertial  acceleration 
vector ,  ^T  ,  is  written  as 

iT  =  a1  x  ^  (b- 

for  a  target  flying  at  a  constant  speed.  Taking  the  cross- 
product  of  the  target's  velocity  vector  into  both  sides  of 
Equation  (B-7)  yields 

Zj.  X  &r=  Vy  X  (li1  X  (B- 

or,  using  the  triple  cross  product  relationship. 


(  B-9  ) 


^  X  ^T=  (V.J.  •  3^)  -  (V^  *  ii'  )  Zp 

Again,  since  the  target's  inertial  velocity  and  angular 
velocity  vectors  are  perpendicular,  the  last  term  on  the 
right-hand  side  is  zero.  Rearranging  Equation  (B-9)  gives 

*  =  (Zp  x  hy)  /  IIV^I2  ( B  - 1 0  ) 

Since 

II tt1  II 2  =  (tt1  *  ii1)  (B-ll) 

substitution  of  Equation  ( B  —  1 0 )  into  Equation  (B-ll)  yields 

ii  li1  ii 2  =  [(JLp  X  AT)*(VT  X  /  HVLpll4  (  B  - 1 2  ) 

In  the  thesis,  then,  the  appropriate  state  estimates  replace 
their  vector  representations  in  Equation  (B-12). 


of  data  for  each  of  the  filters  evaluated  against  each  of  th 
trajectories.  Six  filter/trajectory  combinations  are 
plotted,  in  plot  sections  numbered  one  through  six.  The 
first  three  plot  sections  are  for  each  filter  evaluated 
against  Trajectory  1,  as  defined  in  Section  5.2.  The  last 
three  plot  sections  are  for  Trajectory  2.  Each  plot  section 
is  subdivided  into  plots  a  through  i,  representing  filter 
states  x-^  through  xg,  respectively.  Thus,  the  plotted  data 
is  readily  compared  across  filter  types  and  trajectories  by 
comparing  plots  with  the  same  letter  designator. 


6Cfr  9  828  *  6IJTG  6091  000  0  6091-  6T2X-  628*- 
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The  purpose  of  this  study  was  to  investigate 
target  state  estimation  techniques  for  the  air-to-air 
mode  of  the  AFTI/F-16  automated  maneuvering  and 
attack  system.  The  target  state  estimator  (TSE) 
previously  developed  would  not  perform  to 
specifications;  possible  reasons  for  this  poor 
performance  are  presented  as  well  as  suggestions  to 
upgrade  the  performance. 

This  study  focused  on  three  areas:  (1)  Determine 
if  the  Gauss-Markov  dynamics  model  used  In  the 
current  TSE  was  adequate  for  the  tracking  accuracies 
specified,  (2)  Determine-’' i  f  a  rotation  had  to  be 
performed  to  account  for  the  states  being  expressed 
in  one  frame  while  the  measurements  were  physically 
made  in  another.  ^3)  Determine' what  effect  the 
conventional  covariance  updates,  coupled  with  the 
short  (16-bit)  wordlength  of  the  TSE  computers,  has 
on  the  stability  of  the  Kalman  filter. 

Two  filter  dynamics  models  were  designed, 
tested,  and  compared.  The  first  model  used  complex 
equations  and  closely  modeled  an  air-to-air 
engagement.  Most  of  the  complexity  of  the  model  was 
maintained  in  its  implementation,  and  it  was  used  as 
a  baseline  model.  The  second  filter  used  a  Gauss- 
Markov  dynamics  model  and  made  several  assumptions  to 
simplify  computations. 

Analysis  of  filter  performance  revealed  that  the 
Gauss-Markov  filter  dynamics  model  was,  indeed,  an 
adequate  model.  Also,  the  covariance  matrix  does  not 
have  to  be  rotated  into  the  LOS  frame  if  the 
measurements  are  redefined.  The  second  filter  was 
then  implemented  using  U-D  covariance  factor izat ion 
algorithms,  but  the  time  propagation  routines  used 
were  apparently  flawed.  However,  the  poor 
performance  of  the  TSE  is  no  doubt  caused  by  the 
conventional  Kalman  filter  recursions,  as  they  are 
inherently  unstable. 
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